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Abstract: We consider normal forms in ‘magnetic bottle’ type Hamiltonians of the form H = 
\{p 2 p + ujfp 2 ) + + hot (second frequency u 2 equal to zero in the lowest order). Our main results 

are: i) a novel method to construct the normal form in cases of resonance, and ii) a study of the 
asymptotic behavior of both the non-resonant and the resonant series. We find that, if we truncate 
the normal form series at order r, the series remainder in both constructions decreases with increasing 
r down to a minimum, and then it increases with r. The computed minimum remainder turns to 
be exponentially small in where A E is the mirror oscillation energy, while the optimal order 
scales as an inverse power of A E. We estimate numerically the exponents associated with the optimal 
order and the remainder’s exponential asymptotic behavior. In the resonant case, our novel method 
allows to compute a ‘quasi-integral’ (i.e. truncated formal integral) valid both for each particular 
resonance as well as away from all resonances. We applied these results to a specific magnetic bottle 
Hamiltonian. The non resonant normal form yields theorerical invariant curves on a surface of section 
which fit well the empirical curves away from resonances. On the other hand the resonant normal 
form fits very well both the invariant curves inside the islands of a particular resonance as well as 
the non-resonant invariant curves. Finally, we discuss how normal forms allow to compute a critical 
threshold for the onset of global chaos in the magnetic bottle. 

1 Introduction 

‘Magnetic bottle’ type nonlinear Hamiltonian dynamical systems appear in various areas of 
physics and astronomy. Examples are plasma confinement machines, ion traps and charged 
particle measuring devices, planetary magnetospheres leading, e.g., to the formation of ra¬ 
diation belts, magnetic reconnection, magnetic bottles in the solar corona, etc. (see, for 
example, Dendy 1993, Gurnett and Bhattacharjee 2005, and references therein). 

It is well known that in such configurations there exist both regular and chaotic particle 
orbits. The regular orbits are of oscillatory nature, i.e., a gyration around the magnetic field 
lines combined with a ‘mirror’ oscillation along the field lines. A bouncing of the particles ap¬ 
pears when the magnetic field lines converge towards a preferential direction (see e.g. Jackson 
1962). A magnetic bottle is formed when there are two distinct domains along the field lines 
in which we have such a convergence. The particles are reflected as they approach the two 
‘necks’ of the bottle. The mirror frequency is of order oo z = 0(\V 2 (d 2 Bj_/dpdz)/B z ) 1 /‘ 2 \), 
where V± is the gyration velocity, B z , B± are the measures of the magnetic field along and 
accross the preferential direction (denoted by z) respectively, and p 1 z. One typically has 
w z << w c , where w c = B z q/m is the gyrofrequency of a particle of charge q and mass m. 

A basic form of adiabatic theory (see, for example, Jackson 1962, or Lichtenberg and 
Lieberman 1992), describes mirror oscillations as a consequence of the preservation of a 
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so-called ‘adiabatic invariant’. To lowest order, this corresponds to the particle’s magnetic 
dipole moment p = qVfl / (2 uj c ) . The determination of higher order (in powers of p) adiabatic 
invariants is a classical problem of dynamics. Several methods to deal with this problem are 
discussed in Kruskal (1962), Northrop (1963), Contopoulos (1965), Dragt (1965), Arnold et 
al. (1988), Lichtenberg and Lieberman (1992) and Benettin et al. (1999). 

For magnetic bottle Hamiltonians quite efficient methods can be derived in the context 
of canonical perturbation theory. Such methods have been proposed by Contopoulos and 
Vlahos (1975), and Dragt and Finn (1979). In the canonical context, one starts from the 
basic Hamiltonian 


H= 


a) 


where A is the vector potential corresponding to the magnetic field B = V x A and p are 
generalized momenta conjugate to the position variables. Assuming, in the simplest case, 
axisymmetry around the z-axis, the Hamiltonian in cylindrical co-ordinates (p, z ) takes the 
form 




( 2 ) 


with u i = oj c for particles gyrating around the z-axis. Note that the equations of motion 
arising from the quadratic terms in the Hamiltonian ([2]) represent a case of so-called ‘nilpo- 
tent’ linearization, since one of the eigenvalues of the matrix of the linearized system of 
equations is equal to zero. There is a variety of methods for constructing a normal form for 
such systems (see, for example, Meyer (1984), Cushman and Sanders (1986), Elphick (1988), 
Baider and Sanders (1991), as well as Sanders et al. (2007) for a review). On the other 
hand, the canonical approach leads to a definition of action-angle variables for the magnetic 
bottle, admitting straightforward physical interpretation. In particular, if we define the pair 
of action-angle variables via 


P=\ -sin tpi , p p = \/2wiJi cos , (3) 

V w i 

the quantity J\ is proportional to the magnetic dipole moment, i.e., J\ = mVf[/(2u c ) 
= ( m/q)p . After some steps of perturbation theory, we can then define new canonical 
variables ($i, I±, P^), which are near-identity transformations of the old canonical variables 

((f) i, Ji,z,p z ), so that the Hamiltonian in the new variables takes the form: 

H(e 1 ,h,c,Pc) = z(h,c,Pd + R(Oi,hX,P() ■ ( 4 ) 

The function Z, called normal form, has the form 

Z(h,£,Pc) = +a;|(Ti)C 2 ) + - (5) 

The action I\ would be an exact integral of the Hamiltonian flow under Z alone. On the 
other hand, the function R, called ‘remainder’, depends on the angle Q\, thus it introduces 
some time variations of I\ under the complete Hamiltonian flow of ((5|). Nevertheless, one 
typically has that \R\ « |Z|, implying that the effect of the remainder on dynamics is small. 
Hence, I\ represents a quasi-integral of motion, i.e. a high-order adiabatic invariant, while 
the mirror frequency W 2 (assumed small) is expressed by this theory as a function of I\. 

The convergence behavior of the magnetic bottle normal form series at high normal¬ 
ization orders has not yet been fully explored. A systematic study of the convergence is 
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presented in Engel et al. (1995). These authors considered a polynomial magnetic bottle 
model. Then, they computed a polynomial form of a ‘quasi-integral’ up to degree 14 in the 
canonical variables. From the numerical data, they distinguish three types of behavior, i.e., 
i) convergence, ii) non-convergence, and (iii) ‘pseudo-convergence’, depending on the behav¬ 
ior of the numerical variations of their quasi-integral at various normalization orders up to 
order 14. Here, we extend this analysis to much higher orders, and provide a numerical 
estimate of the asymptotic behavior of the series based on the size of the remainder of the 
normal form. In agreement with basic theory, we observe numerically that the only exist¬ 
ing behavior is ‘pseudo-convergence’, i.e., the series exhibit always an asymptotic behavior. 
This means that an ‘optimal’ normalization order r op t can be identified, up to which the 
remainder decreases in size, yielding the impression that the normalization is convergent, 
while, beyond the optimal order, the size of the remainder increases with the normalization 
order, thus the normalization turns always to be divergent. We also provide evidence that 
the crucial small quantity which enters in all asymptotic estimates is the energy A E of the 
mirror oscillations. Asymptotically, one has the estimates r op t = 0(AE~ a ), for the optimal 
order, and ||i?|| op t = 0(exp(— 1 /AE 7 )) for the size of the remainder at the optimal order, 
with exponents a and 7 specified numerically in section 3 below. We note that theoretical 
exponential estimates on adiabatic invariants in nonlinear modulated oscillators are discussed 
in Neishtadt (1981,1984) and Benettin and Sempio (1994) (for the case of linear oscillators, 
see Howard (1970) and references therein). 

The second main result in the present paper regards the construction of a normal form in 
magnetic bottle Hamiltonians in a case not covered by the usual theory, namely the case of 
resonances. Resonances appear whenever the condition 002/011 = m 2 /mi is satisfied for non¬ 
zero integers mi,m 2 - Such values are marked by the bifurcation of new periodic orbits (in 
pairs stable-unstable) from a so-called ‘central’ (equatorial) periodic orbit (z = 0). The most 
important resonances are of the form uij — nui 2 = 0, with n integer. Resonances of lower and 
lower order appear by increasing the energy. The appearance of the lowest resonances 1:4, 
1:3, 1:2, marks an overall qualitative change of the phase space structure leading eventually 
to the onset of global chaos. 

The usual (non-resonant) normal form can predict the values of I\ when new resonances 
bifurcate, as well as the distance of the resonances from the central orbit as I\ increases 
(Contopoulos and Vlahos 1975). Nevertheless, it cannot describe the structure of the phase 
space near resonances. Here, precisely, we propose a method of construction of a resonant 
normal form for magnetic bottle Hamiltonians, which is applicable both for resonant orbits of 
one (at a time) specific resonance as well as for the non-resonant orbits in its neighborhood. 
Our method can be viewed as a combination of two recently introduced techniques: these 
are i) ‘detuning’ (see Pucacco et al. 2008), and ii) ‘book-keeping’ (see Efthymiopoulos 2008, 
2012). Both techniques reflect ways to optimize the formal treatment of various small quan¬ 
tities appearing in the formal series. The main difference between the usual non-resonant 
construction and the hereby proposed resonant construction is the following: In the case of 
non-resonant series, we select as small parameter either the frequency 002 or the distance (in 
phase space) from the central orbit (see Contopoulos (1965) for a detailed comparison of 
the two approaches). In the resonant series, however, we simultaneously treat the distance 
from the central orbit and the small frequencies as small parameters. We note, finally, that 
the case presently dealt with is quite distinct from the case of resonance in models with a 
periodic space modulation of the magnetic held (as e.g. in Dunnett et al. 1968, McNamara 
1978). 

Implementing our algorithm we compute high order resonant quasi-integrals, and check 
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their degree of accuracy in comparison with the invariant curves found numerically on the 
domain of regular motion, using as reference the same model as Engel et al. (1995). In 
the same way as for the non-resonant normal form, we also here examine the asymptotic 
behavior of the resonant formal series. We demonstrate that in this case as well there hold 
exponential estimates for the dependence of the optimal normalization order, as well as the 
size of the optimal remainder, on the mirror oscillation energy A E. As a final outcome, we 
use the magnetic bottle normal forms in order to analytically compute the critical energy 
beyond which the central periodic orbit becomes unstable. This determines the energy where 
we have the onset of global chaos in the magnetic bottle. 

The paper is organized as follows. Section 2 briefly describes some features of the refer¬ 
ence model (Engel et al. 1995) used in our study, for the paper’s self-containedness. Section 
3 describes the algorithm of computation of the non-resonant normal form. We emphasize 
that this is not a new method but essentially the same as in Dragt and Finn (1979) and 
Engel et al (1995). Also, here as well we employ the method of Lie series which is quite con¬ 
venient for performing near-identity canonical transformations (Hori (1966), Deprit (1969)). 
However, we present the specific algorithmic steps in greater detail than in these earlier pa¬ 
pers, in order to set the context and introduce the notation of the ‘book-keeping’ technique 
(Efthymiopoulos 2012). The main new result in this section regards the numerical study of 
the asymptotic properties of the non-resonant series and the determination of the associated 
exponents entering in exponential estimates of the size of the optimal remainder. These are 
found by extending all normal form computations up to a high order. Section 4 describes the 
novel computation of the resonant normal form construction for magnetic bottle Hamiltoni¬ 
ans. Here, also, we study numerically the normal form’s asymptotic behavior by reaching 
sufficiently high normalization orders. Finally, we compute approximately the threshold for 
the onset of global chaos using normal forms. Section 5 summarizes our conclusions. 


2 Hamiltonian model 


We consider the same axisymmetric magnetic bottle model as in Engel et al. (1995). The 
magnetic held is given by B = V x A, where A is the vector potential given in cylindrical 
coordinates by A = (A pi A^, A z ) = (0, A^,0) with 




( 6 ) 


£>o is the value of the homogeneous (along z) magnetic held component and measures the 
strength of the (octupole) component causing the mirroring effect. In the physical context, 
/3i represents a small parameter. The equations of motion for a particle of mass m = 1 and 
charge q = 1 are derived from the Hamiltonian 


H = - (p — A) 2 = —+ —+ ^ 

2 ; 2 2 2 p 2 p 2 


( 7 ) 


Since </> is ignorable in ([7 pp is an integral of motion. The system can be considered as of 
two degrees of freedom, with ‘effective potential’: 


T j- v\ A (j> (p, Z )pf A\{p,z) 

Veff = V- P — + 


( 8 ) 


Orbits starting with z = 0, p z = z = 0 remain always on the equatorial plane z = 0. For 
such orbits, the equations p = p p = 0 and p p = —dV e ff(p, z = 0)/dp = 0 dehne two types of 
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Figure 1: The family of equipotential curves or curves of zero velocity (CZV), E=V(p,z) = T|/2, 
for various values of E. E cr i t = 0.592 is the value of V(r, 0) at its local maximum, and the dots indicate 
the local minima E = 0. 

equilibria for any fixed radius p: 

i) For Ptf, < 0 the equilibrium solution is p = p c = [—p ( pdA r f > /dp)p =Pc}Z= o] 1 / 2 . This yields a 
circular equatorial orbit surrounding the central axis, with gyration frequency oj c = p,f>/p 2 — 
A<f,(p c ,0)/Pc < 0 . Nearby orbits can be studied by means of the epicyclic approximation. 
Setting i = p-p c and expanding the Hamiltonian around the circular solution we arrive at: 

H = const + ^ + -k 2 £ 2 + (9) 

where k 2 = Sp^/p* + 2 p^A'^/p 2 c - 2p ( pA ( f > /pl - p^A'^/p + A^A^ + ( A (j) ) 2 (derivatives are with 
respect to p, evaluated at p = p c , z = 0). It is easy to check that in Hi the lowest order terms 
quadratic in z are either of order 0{piz 2 ) or 0(^z 2 ), hence small. Thus, the Hamiltonian ([9]) 
is of the general form Q. 

ii) For Pif, > 0 one finds, instead, the equilibrium solution p = p c = p^/A^pc, 0) which 
yields <j) = 0. This solution describes a particle at rest at the distance p c for any value of the 
azimouth. Nearby orbits on the equatorial plane, keeping p^ constant, arise by perturbing the 
radial velocity p ^ 0 while keeping cf> = z = 0. In this case, a similar analysis as above shows 
that nearby orbits describe gyrations around the equilibrium solution, which do not encircle 
the central axis. The gyration frequency is now equal to uj c = k, the minimum and maximum 
distances from the axis are found by the two roots 0 < p\ < p 2 of H = E = V e ff(p, 0), while 
cj) has opposite sign in the intervals p\ < p < p c and p c < p < p 2 - The Hamiltonian is still 
given by ([9]) , and it is easy to check that the lowest order terms quadratic in z are of order 
0(^z 2 ). Thus, the Hamiltonian is again of the general form (|2l). 

The simplest, albeit without loss of generality, case to consider is = 0, which we 
hereafter focus on. In this case, the orbital motion takes place on a meridian plane ( p , z) 
rotating with angular velocity <f> = —A^/p. Gyrating orbits cross the axis p = 0. In order 
to formally avoid discontinuous transitions in the value of 4> at each crossing, the cylindrical 
radius p can be adopted to take both positive and negative values. In units in which Ho = 2, 
Pi = 1 in Eq. ([ 6 ]), the motion on the meridian plane is described by the Hamiltonian: 

H{p, z,p p ,p z ) = ^{p 2 p + p 2 z ) + V(p,z) (10) 
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Figure 2: (a) The main types of orbits for E = 0.2. The central orbit in gray is a typical non-resonant 
regular orbit undergoing mirror oscillations. The ‘figure 8 ’ orbit (black) is a 2:1 resonant orbit. Both 
orbits are limited by the central curves of zero velocity. On the other hand, the left and right gray 
curves show chaotic orbits limited by the outer CZVs. (b) A chaotic orbit for E = 0.7, exploring 
the interior domain limited by the corresponding CZVs. (c) The surface of section ( z,p z ) for p = 0, 
P p > 0 at the energy E = 0.2. Non-resonant regular orbits belong to invariant curves surrounding 
the central (equatorial) periodic orbit at z = p z = 0 , while resonant orbits belong to island chains 
around the origin. The domain of stability is surrounded by a sea of chaotic orbits, (d) The surface 
of section at E = 0.7, where global chaos prevails. 


with a ‘potential’ (equal to A^/ 2 ) given by 


T/v \ 1 2 I 1 2 2 1 4 , 1 2 4 1 4 2, 1 6 

V M=2 P + 2 PZ ~8 P + 8 PZ -W PZ + U8 P 


( 11 ) 


Figure [Tj plots the equipotential curves, or curves of zero velocity (CZV), for different 
values of the energy E = H(p, z) = V(p, z) with p p = p z = 0 (the magnetic field lines have 
a similar form (see figure 1 in Engel et al. (1995)), but form a small angle with the lines of 
Fig[TJ which correspond to constant values of ). These curves define the limits that can 

be reached by the orbits. The potential along the p-axis (z=0) is V(p, 0) = ^p 2 (l — ^-) 2 = E 
and has two minima and one maximum (in each of the half-planes p > 0 or p < 0). The 
maximum is equal to V(p, 0)=E cr jt=16/27« 0.5926 for p= p cr i t =± \/8/3 ~ 1.633. 

Examples of orbits in the above system are given in Fig[2j When 0 < E < E cr a there are 
three permissible regions of orbits, one close to the z-axis (p = 0 ) and two on the right and on 
the left (Fig(2^i). In the central region we find ordered orbits, which obey some quasi-integral 
of motion. The gray and black central orbits in Figj2^i correspond to a non-resonant and 
resonant case respectively. The associated phase portrait (surface of section p = 0, p > 0) is 
shown in Fi^2]:, for the energy E = 0.2. In this case, the main resonances are 2:1 and 3:1, 
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which both produce double chains of islands of stability (four and six respectively in Figj2fc), 
as discussed in detail in section 4. On the other hand, in Figj2^i the left and right orbits in 
gray are chaotic orbits having no intersection with the surface of section of Fig(2jc. Finally, 
for much higher energies (E = 0.7, Fig^ 5 ) 3 ,d), the central orbit has become unstable and the 
system is characterized by global chaos. 

In section 4, we shall employ the normal form method in order to compute the critical 
energy value where the onset of global chaos takes place. 


3 Non-resonant normal form 

3.1 Algorithm 

The Hamiltonian (HOI) is of the general form (jT]) . A non-resonant normal form for this Hamil¬ 
tonian can be constructed by the method of Dragt and Finn (1979) or Engel et al. (1995). 
We summarize here the main steps, introducing our own notation and terminology: 


i) Introduction of complex canonical variables. Introducing the linear canonical change of 
variables 


P = 


q\ + w\ 
\/2wi,o ’ 


Pp 


m + pi 


z = 92, 


Pz = P 2 


( 12 ) 


with cugo equal to the frequency induced by the quadratic term 0 (p 2 ) of the potential (cugo = 
1 in our example), the Hamiltonian takes the form H = H 2 + H± + Hq, with 


#2(91, Pi, P2) = iuiipqiPi + - 


P 2 


(13) 


The terms H 4 and Hq are of fourth and sixth degree respectively. For orbits crossing the 
z-axis (with pp = 0 ), cugo is equal to half the gyration frequency. 


ii) Book-keeping: We organize the terms in the Hamiltonian in groups of ‘different order 
of smallness’. Formally, we introduce a ‘book-keeping’ parameter A, with numerical value 
A = 1, and write the Hamiltonian as 

H = 77 (0) = + A H[ 0) + A 2 H { 2 ] + ... (14) 

The superscript (0) means ‘no normalization step performed so far’. A subscript i, accompa¬ 
nied by a book-keeping coefficient A*, means “i-th order of smallness”. The arrangement of 
all the Hamiltonian terms in the groups Hq*\ h[°\ etc., is done by adopting a ‘book-keeping 
rule’. For polynomial Hamiltonian models, the simplest choice of rule is to associate book¬ 
keeping order with polynomial degree. Thus, in the Hamiltonian (|14D we set to be the 
ensemble of terms of polynomial degree 2s + 2. This renders our non-resonant construction 
equivalent to those of Dragt and Finn (1979) or Engel (1995). However, the above book¬ 
keeping rule will be modified in the resonant construction dealt with in section 4 below. We 
note here that the form of the Hamiltonian (section 2) in the more general case p^ > 0 allows 
for employing the same book-keeping rule as above, with p substituted by the local variable 
£ = p — p c around an equatorial orbit at p = p c . Instead, in the case p§ < 0 we have to add 
one more power of A for every power of the quantity /?i, which now plays also the role of 
small parameter. 
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iii) Choice of ‘kernel set’ A4. To choose a ‘kernel set’ means to answer the question of 
which terms are kept in the normal form along the normalization process. In the present ex¬ 
ample, in the Hamiltonian there appear monomial terms of the form (ifp\ q 2 2 p l 2 - where 
ki, l\, k 2 , 1-2 > 0. For determining a non-resonant formal integral, it suffices to keep 
in the normal form the terms satisfying k\ = l\. Then, the normal form takes the form 
Z = Z(qipi,q 2 ,p 2 ), he. its dependence on qi,pi is only through the product qipi- Then, 
the quantity I\ = iq\p\ is a formal integral, since {Ii,Z} = 0 (where {•,•} denotes the 
Poisson bracket). However, it is practical to exclude some further terms from the normal 
form. Namely, we also exclude the terms with k\ = l\ but I 2 f 2 0, except if k\ = l\ = 0 and 
1 2 = 2, k 2 = 0. This means to retain terms of the form Ifqp, for exponents n / 0, as well 
as the ‘kinetic’ term \p\. Since I\ is a formal integral, after these definitions the normal 
form reduces to the form Z = uj\I\ + p 2 /2 + U(Ii, ( 72 ); he. Z takes the form of an one-degree 
of freedom oscillator (uj\I \) plus a ‘kinetic’ and a ‘potential’ term for the second degree of 
freedom (for which I\ acts as a parameter in the ‘potential’ {/(T), ( 72 ))- 
In summary, as kernel set A4 we choose: 


M = [k] = l] and (k 2 = 0 ,12 = 2 if k\ + l\ = 0, or I 2 = 0 if k\ + l\ > 0)} 


(15) 


iv) Normalization. We construct the non-resonant normal form by using canonical trans¬ 
formations via Lie series (see, e.g., Efthymiopoulos 2012 for a tutorial introduction to Lie se¬ 
ries) . We thus introduce the sequence of transformations: {q\ , (72 , pi , P 2 ) —>• (< 7 ^, q ^, p\ l ^, pP ), 
--X ( q\ \q 2 ,Pi ,p 2 ), —>, ■ ■ •, where superscripts (1), (2),...,denote the new canonical vari¬ 
ables after the first, second, etc., normalization steps. In the r-th step, the transformation 
is given in terms of a Lie generating function Xr { ( lp , qP ■, Pp , pP) ■ The transformation 
of any function f(q^ 1 ^, p 2 l \p\ l \p 2 hi the new variables is found by replacing 

(■qi l \q 2 l \p\ 1 \p2 ^) with (qp, qp,Pi \pP) in the arguments of / and, then, by 
computing f = exp (L Xr )f, where 


OO -1 

eX P( L Xr) = E T\ L Xr 

k =0 


(16) 


L Xr = {-,Xr} denoting the Poisson bracket operator. In particular, the variables themselves 
are transformed according to qp = exp (L Xr )q( E q ^), where q[ r ^(q^) is the identity 
function (and similarly for the remaining variables). In the sequel, for simplicity we drop 
superscripts in the notation of the canonical variables. 

The generating functions y>, r = 1,2,... are computed recursively. Namely, the Hamil¬ 
tonian after r — 1 normalization steps has the form 


H= Z 0 + XZ] + X 2 Z 2 + ... + A r " 1 Z r _i + X r Hp-V + X r+1 HpTp + 


(17) 


where all the terms in Zq, Z\,...,Z r -\ are in normal form, i.e., they belong to X4. According 
to the book-keeping rule choosen in (ii), the terms Z s , or Hs ^ are of polynomial degree 
2 s + 2. In particular, Zq = H 2 = iuii^qipi + p 2 /2. Let Hr denote the terms of Hr 
not belonging to M. The generating function Xr is the solution to the homological equation 

= -A r Hp~V . (18) 


P 2 2 


[Z 0 , Xr} = < iui,oqiPi + y, Ar 


After computing Xr, we compute the transformed Hamiltonian 

= exp (L Xr )H^ . 


(19) 



This is in normal form up to terms of book-keeping order r, namely 


H (r) = Z 0 + AZi + A 2 Z 2 + ... + \ r Z r + A r+1 ^i + A r+2 tf r ( J 2 + ... (20) 


This completes one step of the normalization algorithm. 

In order to solve the homological equation ([18]) . we first decompose Hr- 1 " ^ 
monomials 



E 


hn K 1 n 11 n K2 

n k 1 ,l 1 ,k 2 ,h qi Pl q2 P 2 


feldl ^2 ryJ'2 


ki,li,k2,l2>0 


in the sum of 


k\-\-l\-\-k 2 ~\~l 2 —2 t*+2 


(r—l) 

with known coefficients h ki ^ l . However, one notes that, contrary to the usual Birkhoff 
normal form around elliptic equilibria, in the present case, due to the presence of the term 
p|/2 in Zq, a monomial q'\ l p l \ (/ 2 2 p 2 2 does not constitute an eigenfunction of the linear operator 
D Zo = {Z 0 , •}, since 

Dz^lkq^lk = Kh ~ h)ui,oqi 1 Piq% 2 p l 2 ~ p 1 'p l j +1 . ( 21 ) 


Thus, the homological equation cannot be solved by term-by-term comparison of the coeffi¬ 
cients as in the usual Birkhoff case. Instead, we form the sets of terms 


•A-kl 


(r-l) 

a kl,n 


n k n l a n n 2r+2 ~ k ~ l ~ n 

QlPlQ 2 P 2 


k + l < 2 r + 2, 


n = 0,1,... , 2r + 2 



( 22 ) 


( j *_ (y _ 

where ar kln 1 = —h kln 2 r+ 2 -k-l-n for given values of k = k\ and l = l±. Then, setting the 
generating function \ r to contain a similar group of terms with (yet unknown) coefficients 
b kl rE the homological equation is decomposed in the set of linear systems of equations 


1 c k i 1 

ch 2 



b 

b 


(r—1) 

kl, 0 
h—i) 
kl, 1 


\ ( 


a 

a 


(r~ 1) 
kl, 0 
(r-l) 
kl, 1 


\ 


V 


Cki 2r + 2 
Ckl 


k 


l 




\ 


6 (r_1) / 

u kl, 2 r+ 2 -k-l / 


V 


(r—l) 

a kl, 2 r+ 2 -k-l 




(23) 


where c k i = i(l — k)u\p. If l k, the system (fT~fli can be solved by backward substitution, i.e. 
we first solve the last equation, then substitute and solve the previous one, etc. On the other 
hand, the choice of kernel set M implies that the normalization should eliminate from the 
normal form also terms with k = l and n = 0,1,... 2r+\—k—l. Then we have a kl 2r+ 2 -k-i = 0 
and the last of Eas. (l23l) becomes the identity 0 = 0, while the remaining equations for 

(r—l) 

k = l form a diagonal system with non zero-determinant for the coefficients b kl n ; with 
n = 1,2,... 2r + 2 — k — l in terms of the coefficients a kj ^ with n = 0,1,... 2r + 1 — k — l. 

(r—l) 

Also, the coefficient b kl 0 ' becomes arbitrary, and can be set equal to zero. This completely 
specifies the generating function % r . 
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Figure 3: (a) The Poincare surface of section for E = 0.1 (black curves) superimposed with the curves 
of equal values of the non-resonant formal integral (gray points), truncated at order 5 (corresponding 
to polynomial degree 12). (b) Same as in (a) but for E = 0.2. The non-resonant formal integral 
cannot reproduce the 2:1 islands of stability. 

3.2 Implementation 

Implementing the above procedure to the Hamiltonian (HD, we compute a non-resonant 
normal form using a computer-algebraic program. At low normalization orders, this yields 
results similar to those of Engel et al. (1995). However, we extended the computations up 
to the normalization order r = 15 (corresponding to the polynomial degree 32), while all 
series expressions were truncated at order rt run c = 20 (corresponding to polynomial degree 
42). The asymptotic behavior of the series at high normalization orders is examined in the 
next subsection. At any rate, even at low normalization orders one obtains expressions for 
the non-resonant formal integral comparing well with numerical results for the non-resonant 
ordered orbits. As an example, after five normalization steps, the transformed Hamiltonian 
(restoring the book-keeping constant value A = 1) reads: 

H (5) = h- 0.1875/? - 0.046875/? - 0.0256348/? - 0.0184021/? - 0.0152607/? 

+ 0.5 pj + 0.5hqj + 0.15625 lfq% + 0.1875/? ql + 0.299194/? ql + 0.551285 ifqj 
- 0.151042/? q\ - 0.46224/? q\ - 1.29767/?<?? (24) 

+ 0.107812 l\q\ + 0.669227 /?ql - 0.0697545/?g? + /? (5) 

where I\ = iqipi, and R ^ + H^ + ... denotes the remainder series. The quantity 

/] represents a formal integral of motion, i.e., the non-resonant formal integral. The mirror 
frequency is expressed in terms of I\ by the terms in (|24D quadratic in q 2 - Thus 

wl(h) = h + 0.3125 /? + 0.375/? + 0.598388/? + 1.10257/? + ... (25) 

As explained below, the possibility to compute U 2 in terms of I\ turns to be crucial in the 
subsequent construction of a resonant normal form (section 4). 

In the expression (1241) all symbols refer to the new canonical variables after the com¬ 
position of five canonical transformations, i.e. q\ = q^ etc. Similarly, I\ expresses the 
quasi-integral in the new canonical variables. This can be transformed to the original vari¬ 
ables by computing (up to book-keeping order r) the expression 

$(<?!,Pi, < 12 ,P 2 ) = exp(— L X1 ) O exp(-L X2 )... o exp(-L Xr )(iqipi) . (26) 
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Finally, the integral $ can be expressed in the original variables z,p z by the substitution 

qi = -V2(p-ip p ), pi = -V2{p p - ip), q 2 = z, p 2 = p z • (27) 

This allows to compute theoretical invariant curves on a surface of section p = 0, p p > 0 
for given energy E. To this end, we set p = 0, p p = (2 (E — 17(0,2:)) — p 2 ) 1 / 2 . Then, 
$ is expressed as a function of ( z,p z ) only, i.e. <I> = & sect (z,p z ; E). Figure^ shows a 
comparison between theoretical and numerical invariant curves on the surface of section for 
E = 0.1 (gray points and black curves respectively). The theoretical curves are obtained by 
computing the level curves <& sec t(z,p z ',E) = I ct , where the constant value I c t is computed as 
let = &sect(zo,p z o] E), (zo,p z o) being the initial conditions leading to one invariant curve. In 
Figj3jr the normalization order is relatively low (r = 5). Even so, the theoretical invariant 
curves have a good degree of coincidence with the numerical curves in nearly the whole 
stability domain. In fact, at this energy level no conspicuous resonances are present in 
the interior of the stability domain, while important resonances are only present near its 
boundary. The situation, however, is altered at higher energies, as exemplified in Figj3b, 
referring to the surface of section for E = 0.2. In this case, the numerical surface of section 
exhibits two couples of conspicuous islands corresponding to a (double) 2:1 resonance. As 
shown in section 4, the value of the critical energy E 2: i where the 2:1 periodic orbits bifurcate 
from the center (at z = p z = 0), as well as the existence of a double 2:1 island chain, can 
be predicted already by the non-resonant normal form. We find that the energy E = 0.2 
is a little above the bifurcation value and the 2:1 island chains have moved from the center 
outwards. However, it is obvious that the non-resonant normal form construction is unable 
to reproduce the phase portrait in the zone of the 2:1 resonance. Instead, the non-resonant 
theoretical invariant curves pass through the resonant islands, i.e. they mimic a non-resonant 
behavior. This problem is remedied in section 4 by the construction of a resonant normal 
form. 

3.3 Asymptotic behavior 

The (non) convergence behavior of formal series, in general, is determined by the accumula¬ 
tion, in the series terms at successive orders, of small divisors. In the usual Birkhoff series, 
the pattern of accumulation of divisors can be unravelled by carefully examining the various 
chains of terms produced by the recursive normalization scheme at successive orders (see, for 
example, Efthymiopoulos et al. 2004 for a heuristic analysis of the accumulation of divisors 
in both cases of a non-resonant and resonant Birkhoff normal form around an elliptic equilib¬ 
rium). In the present case of non-resonant normal form, however, the analysis is perplexed by 
the fact that the propagation of divisors depends on the solutions, at each order r, of the non¬ 
diagonal set of Eqs J23p . due to the fact that the original Hamiltonian does not contain a term 
ioj 2 p 2 q 2 . Still, one readily sees that implementing repeatedly Eas. (f23l) at successive orders 
leads to chains of terms growing in size with r by an upper bound Q r (e) = 0(rl a (e b /cof 0 ) r ), 
where a, b, c are positive exponents and e = (|pigi| + \p 2 + q 2 \) 1 ^ 2 is a measure of the distance, 
in phase space, from the origin. This implies that an asymptotic behavior is expected for 
the adiabatic invariant formal series. Hereafter we demonstrate that this is so by numerical 
experiments where, for fixed value of the energy E, we study the behavior of the series at high 
orders as a function of the energy A E of the mirror oscillation, given by A E = E — E\, where 
Ei is the gyration kinetic energy. The energy A E plays here the role of a small parameter, 
since for A E = 0 we have equatorial orbits corresponding to a one-degree-of freedom, i.e. an 
integrable limit of the Hamiltonian model (| 10 1) . 
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Figure 4: Asymptotic behavior of the non-resonant normal form series exemplified by a specific 
computation using the norm definition (1301) for the remainder function with parameters /3 = 0, and 
E = 0.2. (a) Dependence of the truncated remainder norm on the truncation order N at three different 
normalization orders r = 5, r = 8, and r = 13, when A E = 0.001. (b) The quantity ||l?l r ’ 20 )|| as a 
function of r for the indicated values of A E. The optimal remainder corresponds to the minimum 
of each curve, (c) The optimal order as a function of A E in log-log scale. An approximate inverse 
power-law holds for A E below a threshold AE ss 10~ 3 . The fitting curve corresponds to an exponent 
0.15. (d) The value of the optimal remainder (in log | log | scale) versus AE (in log scale). The straight 
line represents an exponential law R opt oc exp((— AEq/AE) d ) with AEq = 10 -3 and d = 0.12. 


12 






To this end, we consider first a norm definition for the remainder series. At the normal¬ 
ization order r, the remainder series has the form (setting A = 1): 

OO OO 

fl w («i,«!2,Pi,P2)= £ Hj">- £ £ • (») 

s=r+l s=r+l fci^i,fc2,^2>0 

kx +/1+A12 +^2=2 s +2 


We define the N-th order truncation of R> r > as 

N N 

R (r ' N) (qi,q- 2 ,Pl,P 2 ) = J2 H s r) = J2 J2 

s=r+l s=r+l ki,h,k 2 ,h >0 

k\ -\-l\-\-k 2 -\-l 2 — 2s-(-2 



(29) 


Let now E, A E be fixed such that 0 < A E < E. Consider a fixed direction q 2 = f3p2 in the 
plane (92,^2)- It is easy to prove that the quantity 




E, AE,P 


s=r+l ki,h,k 2 ,h >0 \ 

k 1 +11 -{- /c2 -f-12—2 s +2 


X 


m 


k 2 


2A E 


&2+ Z 2 

2 


1 + % 2 ((L-AE)M,o)) 


(30) 


satisfies all the properties of norm definition. The norm (|30l) provides a measure of the size 
of the remainder at the given energy levels E, A E. In particular, we find the following: 

i) For A E sufficiently small, the sequence for A = r + 1, r + 2, ... is 

convergent for, N —>• 00, at all normalization orders r. An example is given in FigH^. We 
fix /3 = 0, E = 0.2, A E = 0.001, referring to an estimate of the size of the remainder along 
the axis z = 0 on the surface of section of FigHh, at a distance p z = V2AE « 4.5 x 10~ 2 . 
The behavior of ||i?( r,Ar )|| is shown for three different normalization orders r = 5, r = 8, 
r = 13. In all three cases, we find that the norm of the truncated remainder converges rather 
quickly with the truncation order N. This implies that the value of the remainder found at 
the maximum truncation order N = 20 used here is a good measure of the limiting value 
||^>(r,°o)|| £ or £j iree c h osen normalization orders r. 

ii) Figure [4^ provides an indication of the asymptotic behavior of the series for the 
particular parameters. Namely, we observe that the norm of the remainder decreases as we 
move from the normalization order r = 5 to r = 8, however, it increases as we move from 
r = 8 to r = 13. The dependence of ||i?^ r,Ar ' ) || on r (fixing N to the maximum N = 20), for 
fixed /3, E is shown in detail in FigHJo, for five different values of the mirror oscillation energy 
A E. The asymptotic behavior of the remainder series is evident in this plot, which shows 
also that the optimal order r op t , at which the norm of the remainder becomes minimum, 
decreases as A E increases, while the value of the norm at the optimal order increases with 
A E. 

iii) The dependence of r op t on A E, shown in FigHJz is approximately power-law like. The 
straight line indicates a power law with exponent -0.15 which holds for energies below a 
threshold value A E ~ 10~ 3 . Note that an optimal order as low as r op t = 2 is reached at the 
energy A E ~ 6 x 10~ 2 . This corresponds to p z ~ ±0.35 along the axis z = 0 on the surface 
of section. Simple visual inspection of Fig{2jz shows that this value is close to the border 
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of the stability domain for z = 0. This suggests that the limits of the stability domain can 
be estimated analytically by the requirement that r op t becomes very small., e.g. r opt = 2. 
Physically, this marks the limit of overall validity of the non-resonant normal form. 

iv) The optimal remainder value is exponentially small in 1/A E. Figure 0J1 shows the 
quantity log | log(||i?^ r ° p ‘’ 20 ^||)| vs. log AT 1 (the value of the remainder decreases for higher 
values in the ordinate). The straight line corresponds to a law | |i?( r ° pt )|| ~ exp(— (AEq/AE) d ) 
with d = 0.12. Notice that, as shown also in FigQJx the exponential dependence implies that 
for A E small we can obtain quite small optimal remainder values (e.g. of order 10~ 9 when 
A E = 10 , rising to 10 -4 when A E is of order 10~ 2 . These numbers set the overall level of 

precision of the non-resonant quasi-integrals as a function of the mirror oscillation energy. 

Let us note here that the asymptotic character of the above normal form construction 
is due to the fact that we seek an expansion in which all quantities are defined in an open 
domain of the phase space. If, instead, after a few normalization steps, we fix a value of 
the action I\ = I\*, and expand locally the Hamiltonian (e.g. Eq. (f24l) l around this value, 
we can arrive at a form of the normalized Hamiltonian in action-angle variables, allowing 
for the implementation of the Kolmogorov algorithm of the KAM theorem. The existence of 
invariant curves in the phase portraits indicates that such a process should yield convergent 
series on a Cantor set of initial conditions. However, exploring such convergence is beyond 
the scope of our present study. 


4 Resonant normal form 


4.1 Bifurcation energy 

Resonant periodic orbits m 2 /mi bifurcate from the central equatorial orbit when m 2 /m\ = 
CU2/W1. Using the non-resonant normal form we can predict the bifurcation energy of the 
m 2 -mi family. The gyro-frequency of the equatorial orbit ui :eq is computed by setting q 2 = 
P 2 = 0 in the normal form (as in Eq. (124[) ). Then 


dZ(I 1 ,q 2 =p2=0) 

“ 1 « (/l) = - an - 

The m 2 /mi family bifurcates at the action value /* given by: 

m 2 u)i t e q (Ii) = miuj 2 {l*) 


(31) 


(32) 


where u 2 (h) is computed as in Ea. (f25l) . Finally, the bifurcation energy is computed by 
Em2/ml = Z(It). 


4.2 Algorithm 

Denoting oj* = uq(/*), = U 2 (I*), we now construct a resonant normal form for the m 2 /mi 

resonance as follows: 


i) Book-keeping and Hamiltonian preparation: As long as Ii is considered as a small 
quantity, one has that the quantities wyo — uj\ = O(Ii), uj\ = 0{I\) are also small. Both 
quantities play the role of ‘detuning’ parameters (Pucacco et al. 2008), since they both 
represent a difference from the unperturbed frequencies which are cuyo and zero respectively. 
Fixing a value 7i = /*, we can formally take this fact into account in the Hamiltonian by 
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adding and substracting the above quantities with a different book-keeping factor , i.e. we set: 

H = iwxxi'Pm + t 2 pI + ■■■ = iu*Piqi + ^pi + -(w| )V 

—A ^i(u;* — + -iaifff'z 2 + H^j + 0(A 2 ) (33) 

Since A = 1, nothing has really changed with respect to the original Hamiltonian. How¬ 
ever, the second frequency was now explicitly introduced in the zero order Hamiltonian term 
which is subsequently used in the normal form construction. Furthermore, if cuf and iwf 
satisfy a resonant condition, it is possible to proceed with the resonant form of the Birkhoff 
normal form. We note that the correspondence between book-keeping orders and polynomial 
degrees is now broken, namely, at the book-keeping order r one has terms of the polynomial 
degrees 2,4,..., 2r + 2. However, this poses no formal obstacles to the construction of the 
normal form. 0 

The remaining steps in the normal form construction are standard (see Efthymiopoulos 
2012 for a review). Introducing the linear canonical change of variables 


Q2+W2 \/wI(h?2 + P2) 

~~ y/%% ’ PZ ~ V2 

the Hamiltonian becomes of the general form (1141) with 

Hq 0) = iulqipi + • 

The functions contain monomials of the form q^'pfq^pf- 

2 < Aq + &2 T /i -M 2 < 2 ?’ + 2 . 

ii) Choice of kernel set M. For the resonance m 2 /mi we set 


(34) 

(35) 

The book-keeping rule is 


■Mres = {(Aq, k 2 ,h, I 2 ) such that (Aq — li)m\ + (Ai 2 — h)^a 2 = 0} . (36) 


This choice ensures that the quantity I res = i{rn\qip\ + ni2q2P2 ) is a formal integral in the 
new canonical variables. 

iii) Normalization. The normalization proceeds recursively by the same scheme as in 
subsection 4.1. The functions Hr 1 denote now the terms of ' not belonging to A4 res . 
The homological equation at the r-th step has the form 

{iu{qiPi + w\q 2 P2, Xres,r} = -H/f^ . (37) 

Thus, the equation is diagonal with respect to all monomials belonging to Hf , i.e., for 
every monomial term ; 2 if 1 p l \ (-Pip/ i n Hr we add a corresponding term in Xres,r 

with coefficient h^~j ^ k2 l 2 /(i((ki - h)u* + (k 2 - h)^)). 

An alternative algorithm to the above would be to introduce a local action variable 
J\ = I\ — If to be treated as a new small parameter. However, in practice we found that this 
approach has worse convergence properties than the ones induced by the normal form after 
the manipulation of the Hamiltonian as in Eg. (1331) . The latter’s efficiency can be tested by 
numerical examples as below. 

1 The following comment offers some insight into the whole above ‘book-keeping’ process: in the case of the 
2:1 resonance treated in detail below, we find |(coj —wi,o)| ~ 0.1, (oij) ~ 0.25. The first quantity can be called 

“of order 1”, but, for the second, the characterization “order 0” or “order 1” would be equally acceptable in 
practice. Note, however, that the fourth order terms in the Hamiltonian include a term h. 22 p 2 z 2 , with a real 
coefficient /122 and p ~ (2Ji) 1//2 . One then finds that the value of /122 has to be such that at the particular 
radius p*, corresponding to the resonant value /*, one will obtain that /?22 (p* ) 2 has only a small difference 
(‘of order one’) from 2(u4 j) 2 )- This is reflected by our choice of book-keeping, which results in the quantity 
NJ 122 P 2 — {1/2){lo2) 2 ]z 2 formally appearing in the Hamiltonian (1331) . 
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Figure 5: The Poincare surface of section ( z , p Z: p = 0) (black curves) superimposed with the curves 
of equal values of the integral of the resonant normal form (gray points) for (a) E = 0.115 and (b) 
E = 0.20. The 3 : 1 resonant islands (two sets of 3 islands) in (a) and the 2 : 1 resonant islands (two 
sets of 2 islands) in (b) are very well reproduced by the resonant normal form. 


4.3 Implementation 

The resonant formal integral I re s can be transformed to an expression in the original variables 
in the same way as for the non-resonant case. After r normalization steps we have 

$res(qi,Pi,<i2,P2) = exp(-L Xres ,) o exp(-L Xres2 )... o exp(— L Xresr )(irri\q\p\ + im 2 q 2 p 2 ) ■ 

..... (38) 
This allows, again, to compare theoretical with numerical curves on the Poincare surface of 

section. Figure [5] shows an example, referring to two different resonances, namely 3:1 (FigJSfe.) 

at the energy E = 0.115, and 2:1 (FigJSJj), at the energy E = 0 . 20 , same as in Figj3b- In 

both cases, the energy is taken close to but above the bifurcation energy value, which, using 

the normal form approach (see subsection 4.1) at the 8 -th order, is found to be -E 1/3 = 

0.097279 and E x / 2 = 0.188036 (the values determined numerically are Ei/ 3 = 0.097253 and 

E \/2 = 0.188015 respectively). In both cases the resonant normal form represents well the 

corresponding islands of stability. Furthermore, in both cases the resonant normal form 

represents also well the invariant curves which surround the center both in the interior and 

the exterior of the resonant zone. The accuracy of the normal form computations, which 

depends on the behavior of the remainder of the normal form series as a function of the 

normalization order, is examined in detail in the next subsection. Here, we emphasize that 

the overall limit of validity of the resonant normal form approach is defined by the appearance 

of other resonances, of different order than the resonance under consideration. These extra 

resonances are conspicuous in the outer parts of the surface of section of FigsO^gb, where we 

see also the beginning of a chaotic sea surrounding the main domain of stability. 

4.4 Asymptotic behavior 

The asymptotic behavior of the resonant normal form is probed again by numerical experi¬ 
ments, in the same way as in subsection 4.3 for the non-resonant normal form. Here, since 
both frequencies uif and uj 2 * are fixed, we introduce a slight modification of the norm defi¬ 
nition with respect to Eg. (13011 . taking into account also the different choice of book-keeping 
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Figure 6: Asymptotic behavior of the resonant normal form, (a-d) Same as in Fig]4jr-d, but for the 
resonant normal form computation with parameters as indicated in the figure. The norm definition 
is given in Eci. (l39l) . 


rule. Thus, we set: 


II * 0 


r,N) | 


N 

e,ae,p = E 

s=r-\-l 


E 

ki,li,k2,h>0 

&i+^i+fc2+^2=2,4,...2s+2 




k -+■ /1 


I&2 


2A E 


1 + P 2 (uj%) 2 ((E - AE)/ut)) 


k 2 + l 2 

2 


(39) 


Figure [Gj which is quite similar to FigdJ clearly shows that the resonant normal form series 
exhibit asymptotic properties analogous to the non-resonant one. Nevertheless, a comparison 
of Figs. and [6}o shows that the overall error of the resonant formal integrals is uplifted 
by about one order of magnitude with respect to the non-resonant case for similar levels of 
mirror oscillation energy A E. The exponent found in Fig. 8c is also close to 0.2. Finally, 
we observe that the exponential regime for the scaling of the optimal remainder with 1/A E 
holds for mirror oscillation energies below a value A E < 10 -3 . 


4.5 Normal form determination of the threshold to global chaos 

The domain of regular motion around the central periodic orbit shrinks, in general, as the 
energy increases. This domain shrinks to zero at an energy Et where the central orbit 
exhibits (for the first time) a transition from stability to instability. The energy Et can 
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Figure 7: (a) Theoretical values for the frequencies lo\ and W 2 as functions of the energy E by the 
normal form at the normalization order r = 10. The intersection point yields an estimate of the 
energy E t where the first transition of the central periodic orbit from stability to instability takes 
place, (b) The theoretical estimate E t as a function of the normalization order r. The horizontal line 
corresponds to the numerically computed value of E t . 


be found numerically by computing the monodromy matrix of the central periodic orbit at 
different energies. Numerically, we find Et = E\ ~0.36688. 

The energy Et marks the limit of applicability of the normal forms in magnetic bottle 
problem, in either the non-resonant or resonant form. We show now how the energy Et 
can be determined by normal form computations. We proceed as follows: at the energy 
E = Et, the two frequencies uj± and u >2 become equal (Contopoulos, 1968). We then have the 
bifurcation of two equal period resonant periodic orbits from the central orbit. Thus we have 
Et = E]_j\. However, E^n can be computed as indicated in subsection 4.1, for the resonance 
m\ = m 2 = 1. 

The accuracy of the theoretically computed value depends on the normalization order r. 
FigEh shows the frequencies u\ and C 02 as functions of the energy E using the normal form at 
the normalization order r = 10. The intersection point of the two curves yields a theoretical 
estimate £)=0.39550, which has an error 5E ~0.0286. The error is reduced as r increases. 
Figure [TJo shows the estimate for Et as a function of r up to r = 30. The convergence to the 
numerically computed value is rather slow, the error being about 10~ 2 at r = 30. 

5 Conclusions 

In the present paper we explored the limits of applicability of the normal form theory in 
a polynomial magnetic bottle Hamiltonian model. We focused on a new algorithm for the 
construction of the normal form and the computation of quasi-integrals (truncated formal 
integrals) in cases of resonance between the gyration and the mirror frequencies. Furthermore, 
we explored the asymptotic behavior of both the non-resonant and the resonant normal form 
series. Our main conclusions can be summarized as follows: 

1) We explored the asymptotic behavior of the non-resonant normal form series. Extend¬ 
ing the computations at high normalization series confirms the basic theoretical picture that 
the behavior of the normalization is asymptotic. Namely, although the size of the remainder 
of the normal form series goes to infinity when the normalization order r tends to infinity, 
we observe that initially (at low order r) the size of the remainder decreases with r. The 
remainder becomes minimum at an optimal order r which scales approximately as an inverse 
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power of the mirror oscillation energy A E. The size of the optimal remainder is found to 
be exponentially small in 1/A E. We estimate numerically the exponents related to this 
asymptotic behavior. 

2) The non-resonant normal form allows to compute theoretically the energies at which 
resonant periodic orbits of any resonance m 2 : rri\ (between the mirror and the gyration 
frequencies, with mi, m 2 integers) bifurcate from a ‘central’ (equatorial) orbit. We propose 
a novel computation of the normal form in the case of resonances. This is based on combin¬ 
ing two algorithmic techniques called ‘detuning’ (Pucacco et al. 2008) and ‘book-keeping’ 
(Efthymiopoulos 2012). We give numerical examples of applicability of the resonant formal 
series in the case of the 3:1 and 2:1 resonances, and demonstrate their ability to predict the 
form of the phase portrait in the neighborhood of each resonance. 

3) We explore the asymptotic behavior also of the resonant formal series, which is found 
to be qualitatively similar to the non-resonant case, and estimate numerically the associated 
exponents. 

4) The suggested normal form computations serve to predict two results regarding the 
onset of chaos, i) At low energies, one can estimate the limits of the domain of stability 
around the central equatorial orbit, i.e. how far from this orbit (in phase space) does chaos 
become important, ii) The onset of global chaos can be approximated by the energy value 
where the central orbit suffers its first transition from stability to instability, which coincides 
with the bifurcation energy of the 1:1 resonance. 
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